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Abstract 

Current metabarcoding studies aiming to characterize microbial communities generally rely on the amplification and 
sequencing of relatively short DNA regions. For fungi, the internal transcribed spacer (ITS) region in the ribosomal RNA 
(rRNA) operon has been accepted as the formal fungal barcode. Despite an increasing number of fungal metabarcoding 
studies, the amplification efficiency of primers is generally not tested prior to their application in metabarcoding studies. 
Some of the challenges that metabarcoding primers should overcome efficiently are the amplification of target DNA strands 
in samples rich in non-target DNA and environmental pollutants, such as humic acids, that may have been co-extracted with 
DNA. In the current study, three selected primer pairs were tested for their suitability as fungal metabarcoding primers. The 
selected primer pairs include two primer pairs that have been frequently used in fungal metabarcoding studies (ITS1F/ITS2 
and ITS3/ITS4) and a primer pair (ITS86F/ITS4) that has been shown to efficiently amplify the ITS2 region of a broad range of 
fungal taxa in environmental soil samples. The selected primer pairs were evaluated in a 454 amplicon pyrosequencing 
experiment, real-time PCR (qPCR) experiments and in silico analyses. Results indicate that experimental evaluation of 
primers provides valuable information that could aid in the selection of suitable primers for fungal metabarcoding studies. 
Furthermore, we show that the ITS86F/ITS4 primer pair outperforms other primer pairs tested in terms of in silico primer 
efficiency, PCR efficiency, coverage, number of reads and number of species-level operational taxonomic units (OTUs) 
obtained. These traits push the ITS86F/ITS4 primer pair forward as highly suitable for studying fungal diversity and 
community structures using DNA metabarcoding. 



Citation: Op De Beeck M, Lievens B, Busschaert P, Declerck S, Vangronsveld J, et al. (2014) Comparison and Validation of Some ITS Primer Pairs Useful for Fungal 
Metabarcoding Studies. PLoS ONE 9(6): e97629. doi:10.1371/journal.pone.0097629 

Editor: Brett Neilan, University of New South Wales, Australia 

Received January 17, 2014; Accepted April 21, 2014; Published June 16, 2014 

Copyright: © 2014 Op De Beeck et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits 
unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 

Funding: The study was made possible by funding of the Research Foundation Flanders (FWO). The funders had no role in study design, data collection and 
analysis, decision to publish, or preparation of the manuscript. 

Competing Interests: The authors have declared that no competing interests exist. 
* E-mail: jan.colpaert@uhasselt.be 



Introduction 

Until the late 1980s, microbial ecologists and taxonomists relied 
on culturing and morphological and physiological characteristics 
to describe microbial communities and members thereof. In the 
last two decades, DNA sequencing has revolutionized the way 
microbial communities are being characterized [1,2]. In addition, 
since the introduction of pyrosequencing by Margulies et al. [3] , 
characterization of microbial communities has undergone a 
second revolution as this technology (used by e.g. Sogin et al. 
[4] and Buee et al. [5]) enables detailed microbial community 
characterization at greater sequencing depth than was deemed 
possible via cloning and Sanger sequencing. A number of next- 
generation sequencing technologies now enable researchers to 
identify a large number of organisms from environmental samples 
using relatively short DNA sequences. This molecular identifica- 
tion method has been termed metabarcoding [6]. Nevertheless, 
whatever sequencing technology is used, DNA metabarcoding 
generally depends on the amplification of barcode regions using 
taxon-specific primers [7]. Such primers need to be universal 
enough to cover a large group of taxa (e.g. the fungal kingdom), 
but at the same time have to result in amplicons that are variable 
enough to efliciendy distinguish between closely related species or 



to identify operational taxonomic units (OTUs) [7,8]. For fungi 
and oomycetes, the internal transcribed spacer region (ITS; 
spanning the ITS1, 5.8S and ITS2 regions) in the ribosomal RNA 
(rRNA) operon has been recognized as the formal DNA barcoding 
region [9-1 1]. 

The full ITS region in fungi has an average length of 500 and 
600 base pairs (bp) for ascomycetes and basidiomycetes, respec- 
tively, and an average length of 600 bp across all fungal lineages 
[12]. As current 454 amplicon pyrosequencing (using Roche's 
Genome Sequencer FLX (GS-FLX) instrument and Titanium 
chemistry) generates read lengths averaging 450 bp, it is impos- 
sible to span the entire ITS region in a single run. Even with recent 
advances in sequencing technologies that enable sequencing across 
the entire ITS region, it will probably remain desirable for fungal 
metabarcoding studies to exclude the 5.8S region of the rRNA 
operon. The inclusion of conserved regions in DNA sequences is 
known to increase the risk of chimera formation during PCR [13]. 
Therefore, generally, either the ITS1 or the ITS2 region is used in 
ecological studies aiming at the characterization of fungal 
communities. 

Primers that will be used in metabarcoding studies should be 
able to efficiently amplify their target DNA regions in the presence 
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of high concentrations of non-target DNA and contaminants, such 
as humic acids, that may have been co-extracted with DNA [14]. 
Therefore, in silico testing of primers is expected to result in an 
incomplete picture of how primers will behave during amplifica- 
tion of DNA extracted from environmental samples. Comparing 
the amplification efficiency and robustness of primers used in 
metabarcoding studies is important because differences in primer 
efficiency may result in strong biases in favour of more easily 
amplifiable sequences during PCR reactions, potentially influenc- 
ing our view on fungal communities [15—17]. Moreover, a primer 
set that covers a large proportion of the species that compose a 
community of interest and that produces a reliable outcome is 
desired as ecological metabarcoding studies typically rely on a 
single primer pair to map microbial diversity. 

The most commonly used primers in fungal ecology for 
sequence-based fungal identification at the species level were 
published by White et al. [18]: ITS1, ITS2, ITS3 and ITS4, and 
by Gardes and Bruns [19]: ITS IF and ITS4B. Whereas the 
primers developed by White et al. [18] had a broad spectrum, 
ITS IF and ITS4B were developed to be specific for fungi and 
basidiomycetes respectively [19]. ITS IF is most frequendy 
combined with ITS2 to amplify the ITS1 region of the fungal 
rRNA operon and ITS 3 is usually combined with ITS4 to amplify 
the ITS2 region. These primer pairs have been used in many 
branches of mycological research in the past twenty years and are 
popular tools in recent fungal community studies as well [5,20-24] 
(also reviewed in Hibbett et al. [25]). 

The aim of the current study was to evaluate the amplification 
efficiency of these established primer pairs and to compare them to 
a selected primer pair (ITS86F/ITS4) that has been shown to 
specifically and efficiendy amplify ITS sequences from a broad 
range of fungal taxa in human blood samples as well as in 
environmental soil samples [26]. 

Materials and methods 

Study site and soil sampling 

A pioneer pine forest on a stabilised sand dune in the northern 
part of Limburg, Belgium (Hechtel-Eksel: 51°7'33"N, 5°22'22"E) 
was selected to obtain samples for this study. The study site is not 
freely accessible. To gain access to this study site, please contact 
the responsible authorities (Table SI). The soil in this study site is a 
dry sandy soil without a litter layer, poor in organic matter and 
slightly acidic. The average organic carbon content for this site is 
0.7% and the average pH is 4.7. The pioneer vegetation at the 
study site is dominated by young Scots Pine trees {Pirns sybestris L.), 
mosses and lichens, with only few grasses and heather shrubs 
(Calluna vulgaris (L.) Hull). Tree ages at the time of sampling ranged 
from one to five years. The region has an average annual rainfall 
of 800 mm per square meter and the average annual temperature 
is 10°C (Royal Meteorological Institute, Ukkel, Belgium). 

Soil samples for fungal community characterization were 
collected in November 2009. Samples were collected at a depth 
of 0 to 20 cm using a soil corer with a diameter of 1 cm. Four 
replicate soil samples were collected within a distance of ten 
centimetres from each other for seven sampling locations. Each 
sampling location was chosen close to a three to five year old pine 
tree randomly selected in the field. Selected pine trees were at least 
20 m apart from each other. The 4 replicate soil samples were 
pooled for each sampling location, resulting in a total of seven 
pooled samples. Samples were sealed in plastic bags and tighdy 
closed to prevent desiccation during transportation. Upon arrival 
in the lab, soil samples were sieved using a 2 mm sieve to 
homogenize the sample and remove roots, large pieces of organic 



matter and stones. Samples were subsequently stored at — 80°C 
until DNA was extracted. No protected species were sampled 
during the study. 

DNA extraction, PCR amplification and pyrosequencing 

Approximately 250 mg of soil was used for each DNA 
extraction. DNA was extracted in quadruplicate from each pooled 
sample using the UltraClean Soil DNA Isolation Kit according to 
the manufacturer's protocol (MoBio, Carlsbad, CA, USA). This 
resulted in four replicates for each of seven pooled soil samples. 
Subsequently, amplicon libraries were created using barcode- 
tagged primers for the primer pairs ITS1F/ITS2, ITS3/ITS4 and 
ITS86F/ITS4 (Table 1). Both forward and reverse primers were 
synthesized with a tail containing the Roche 454 pyrosequencing 
adaptors and a sample-specific 10 bp barcode (multiplex identi- 
fiers: MIDs) [28] enabling sorting out the obtained sequences after 
sequencing (Roche Applied Science, Mannheim, Germany). 
Fusion primers were designed according to the scheme provided 
in Table S2. 

DNA samples were amplified using a Techne TC-5000 
thermocycler (Bibby Scientific Limited, Staffordshire, UK) under 
the following conditions: initial denaturation at 95°C for 
2 minutes, followed by 40 cycles of denaturation at 95°C for 
30 seconds, annealing at 55°C for 30 seconds and extension at 
72°C during 1 minute; a final extension phase was performed at 
72°C during 10 minutes. Reactions were carried out in 25 ul 
reaction volumes using the FastStart High Fidelity PCR System 
(Roche Applied Science, Mannheim, Germany). Each reaction 
contained 2.75 ul FastStart lOx reaction buffer, 1.8 mM MgCl, 
0.2 mM dNTP mix, 0.4 uM of each primer, 1.25 U FastStart 
HiFi polymerase and 5 ng template DNA (as measured by a 
Nanodrop spectrophotometer). 

Amplified DNA was cleared from PCR primers and primer 
dimers using the Agencourt AMpure XP System according to the 
manufacturer's protocol (Beckman Coulter, Brea, CA, USA). 
Finally, purified dsDNA was quantified with the Quant-iT 
PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, CA, USA) 
and a Fluostar Omega plate reader (BMG Labtech, Ortenberg, 
Germany) and subsequendy pooled in equimolar concentrations. 
The resulting amplicon pool, containing all 84 samples, was 
sequenced on one fourth of a Pico Titer Plate on a Roche Genome 
Sequencer FLX System using Titanium chemistry (Roche Applied 
Science, Mannheim, Germany) according to the manufacturer's 
instructions. 

Bioinformatics processing 

The standard flowgram format (SFF) file that resulted from the 
interpreted flowgrams was deposited in the NCBI Sequence Read 
Archive under accession number SRP026207 (SRA, http:/ /www. 
ncbi.nlm.nih.gov/Traces/sra). From the original SFF file, three 
separate quality and fasta files were created with a custom Python 
script according to the three primer pairs used (Table S2). Further 
analyses were carried out in Mothur 1.31.2 on the individual fastq 
and fasta files [29]. Quality trimming in Mothur removed reads 
shorter than 200 bases, reads longer than 600 bases, reads with 
homopolymers longer than 8 bases and reads containing 
ambiguous bases. Reads were trimmed when the average Phred 
quality score dropped below 35 over a window of 50 bases. Next, 
sequences were compared to each other and duplicate sequences 
were replaced by a single sequence, while archiving the abundance 
data of the unique sequences. Subsequently, unique reads were 
checked for chimeric sequences using the Uchime tool in Mothur 
followed by their removal from the datasets. Unique reads were 
aligned with the pairwise alignment tool in Mothur. Finally, 
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Table 1. Primers used in the current study. 





Primer name 


Primer sequence (5 -3 ) 


rRNA operon binding site 


Reference 


ITS IF (F) 


CTTGGTCATTTAGAGGAAGTAA 


Small subunit 


[19] 


ITS2 (R) 


GCTGCGTTCTTCATCGATGC 


5.8S 


[18] 


ITS3 (F) 


GCATCGATGAAGAACGCAGC 


5.8S 


[18] 


ITS4 (R) 


TCCTCCGCTTATTGATATGC 


Large subunit 


[18] 


ITS86F (F) 


GTGAATCATCGAATCTTTGAA 


5.8S 


[27] 


ITS86R (R) 


TTCAAAGATTCGATGATTCAG 


5.8S 


[26] 



Primers are indicated as forward (F) or reverse (R). ITS86R contains a wrong base at the 3' end. The G should be replaced by a C (see Discussion). 
doi:1 0.1 371 /journal.pone.0097629.t001 - 



species-level OTUs were defined based on a 97% sequence 
similarity level, which is within the range of intraspecific ITS 
sequence similarity [30]. In order to further remove potential 
sequencing errors from the analysis, global singletons (i.e. OTUs 
represented by only a single sequence over an entire dataset) were 
removed [24]. 

Because the primer pairs resulted in different amounts of reads 
per sample, the number of reads per sample were rarefied to 200 
reads per sample. Samples for which less than 200 reads were 
obtained were removed from the dataset. For ITS1F/ITS2 14 of 
28 samples were removed. For ITS3/ITS4 4 samples were 
removed and for ITS86F/ITS4 no samples were removed. Inter- 
sample rarefaction curves were constructed based on 10,000 
iterations. Subsequently, intra-sample diversity, richness and 
Good's coverage estimates were calculated in Mothur 1.31.2 
based on 10,000 iterations. BLAST searches for a representative 
sequence of each OTU (as determined by Mothur) were 
conducted using the PlutoF v2.0 massBLASTer online tool [31]. 
Reads were blasted against the UNITE [32] and INSD [33] 
databases. Resulting HTML files were combined with the 
abundance data obtained in Mothur using a custom Python 
script. This script also acquired the names of species or genera that 
resemble Latin binomials with the highest BLAST score, avoiding 
unidentified OTUs in the databases to be seen as best BLAST hits. 
Unidentified OTUs were indicated as "not applicable (NA)". 

Quantitative real-time PCR 

To evaluate the performance of the primer pairs amplifying 
target DNA from a heterogeneous pool of DNA in environmental 
samples, all primer pairs were tested in a qPCR set-up. A 2-fold 
dilution series (1:1 to 1:64) was made from twelve DNA samples 
(ranging from 5 ng J_ll 1 to 78 pg ill , including one no-template 
control (NTC) for each sample). Amplification was performed in 
optical 96-well plates using a 7500 Fast Real-Time PCR System 
(Applied Biosystems, Foster City, CA, USA) and SYBR Green 
chemistry. PCR conditions were as follows: initial denaturation at 
95°C for two minutes, followed by 40 cycles of 95°C (30 s), 55°C 
(30 s) and 72°C (60 s) and a final extension phase at 72°C for 
1 0 minutes followed by the generation of a dissociation curve to 
verify amplification specificity. These qPCR conditions were 
chosen to mimic the PCR conditions used during the PCR step 
prior to emPCR and amplicon pyrosequencing. Reactions 
contained 2.5 |lL template DNA, 5 uL 2 x Fast SYBR Green 
Master Mix (Applied Biosystems, Foster City, CA, USA), 0.3 ul 
forward and reverse primers (3.3 |lM each) and 1.9 |lL nuclease- 
free H z O in a total volume of 10 (iL. PCR efficiencies (E) were 
calculated as E = (10" 1Alc>pt - 1) x 100. 

To assess a potential PCR-bias at the phylum level, DNA was 
extracted from 15 pure cultures provided by the Mycotheque de 



FUniversite Catholique de Louvain (BCCM/MUCL) including 5 
basidiomycetes (Lentinula edodes (MUCL 44827), Agrocybe praecox 
(MUCL 46727), Coniophora marmorata (MUCL 39471), Suillus luteus 
(UH-Slu-LM8-nl) and Antrodia vatikntii (MUCL 54533)), 5 
ascomycetes [Cladosponum cladosporioid.es (MUCL 53652), Cryptospor- 
iopsu radicicola (MUCL 53485), Monilinia laxa (MUCL 30841), 
Arthroderma otae (MUCL 39756) and Galactomyces geotrichum (MUCL 
52377)), 2 glomeromycetes (Rhizophagus clareus (MUCL 46238) and 
Rhizophagus sp. (MUCL 41833)) and 3 zygomycetes (Mortierella 
verticillata (MUCL 9658), Absidia corymbifera (MUCL 38907) and 
Mucor hiemalis (MUCL 15439), also see Table S4). DNA was 
extracted from cultures using the DNeasy Plant Mini Kit 
according to the manufacturer's instructions (Qiagen, Venlo, 
Netherlands). DNA concentrations extracted from pure cultures 
used for qPCR ranged from 5 ng ju.1 - 1 to 20 ng u\l~ 1 . PCR bias at 
the phylum level was tested according to the qPCR protocol 
described above. 

In silico evaluation of primer pairs 

To evaluate the primer-to-target mismatches in silico, primers 
were tested with PrimerProspector 1.0.1 [34] against sequences 
downloaded from NCBI. Three sets of sequences were download- 
ed from NCBI containing only full-length fungal 5.8S, 18S and 
28S sequences. Duplicate sequences were removed using Mothur 
1.31.2. ITS1F was tested against 3,748 18S rDNA sequences. 
ITS2, ITS3 and ITS86F were tested against 4,421 5.8S rDNA 
sequences. ITS4 was tested against 4,270 28S rDNA sequences. 
For comparison, also all primers investigated by Ihrmark et al. 
[35] and Toju et al. [36] were tested [18-19,27,35-38]. All tests 
were performed as described by Walters et al. [34] using standard 
settings. Primer scores were calculated based on the following 
formula: weighted score = non-3' mismatches xO.40+3' mis- 
matches xl.00+non-3' gaps x 1.00+3' gaps x3.00. An additional 
penalty score of 3.00 was assigned if the final 3' base of a primer 
had a mismatch with its target sequence [34] . 

Statistical analysis 

Statistical analyses were conducted in R 2.13.0 (The R 
Foundation for Statistical Computing, Vienna, Austria). Normal 
distributions of the residuals of models were checked with the 
Shapiro-Wilk test, while homoscedasticity of variances was 
analysed using either Bartlett's or the Fligner-Killeen test. 
Depending on the distribution of the estimated parameters, either 
ANO VA or the Kruskal- W allis Rank Sum Test was used to check 
for significant differences in variances of parameters. Two-by-two 
comparisons were conducted using either Tukey's Honest 
Significant Differences tests or Pairwise Wilcoxon Rank Sum 
Tests. Poisson corrections were implemented for abundance data. 
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Distributions of ratios were compared with Pearson's Chi-squared 
tests. Non-metric multi-dimensional scaling (NMDS) was per- 
formed using the Vegan 2.0 - 8 package in R. 

Results 

Parametrical analysis of 454 amplicon pyrosequencing 
data 

For the three tested primer pairs, GS-FLX sequencing of the 
amplicon libraries generated a total of 151,650 reads. For a read to 
be successfully assigned to a sample, we required that both the 
forward and the reverse MIDs and primers were identified in a 
read with no more than one erroneous base in the MIDs and no 
more than two erroneous bases in the primer sequences. Based on 
the primer and MID sequences, 65,133 reads were assigned to 
their respective sample and 86,517 reads remained unassigned. 
The average length of reads assigned to either ITS1F/ITS2, 
ITS3/ITS4 or ITS86F/ITS4 prior to quality checking and 
trimming was 314, 331 and 369 bp respectively (excluding 
primers). The average read length of the unassigned reads was 
1 16 bp (including primers, data not shown). 

Rarefaction curves were constructed showing the rarefied 
number of OTUs defined at a 97% sequence similarity threshold 
relative to the number of samples (Fig. 1). These results indicate 
that, on average, a higher OTU richness and a better coverage of 
the fungal community can be expected for the ITS86F/ITS4 and 
ITS3/ITS4 primer pairs. The lowest OTU richness and coverage 
was predicted for the ITS1F/ITS2 primer pair. As most 
rarefaction curves tended towards saturation, the sequencing 
depth was assumed to be sufficient to retrieve the most abundant 
fungal OTUs in analysed soil samples that are detectable by the 
respective primers and 454 amplicon pyrosequencing. 

To compare primer pair performance in the 454 amplicon 
pyrosequencing experiment, averages of the number of reads were 
calculated across replicates (four replicates per sample) and 
samples (seven samples) for each primer pair. The average 
number of reads per sample obtained by ITS1F/ITS2, ITS3/ 
ITS4 and ITS86F/ITS4 after quality trimming differed signifi- 
cantly (p<0.01) and primer pairs yielded on average (± standard 
error) 356 (±26), 523 (±43) and 797 (±34) high quality reads per 
sample, respectively (Fig. 2A). The average number of OTUs 
found for each primer pair at a 97% sequence similarity threshold 
(observed OTU richness) also differed significantly (p<0.01). The 



400 i 




Number of samples 

Figure 1. Rarefaction curves for each of the three primer pairs used in this study: ITS1F/ITS2, ITS3/ITS4 and ITS86F/ITS4. In these 
graphs, the number of samples is plotted against the rarefied number of operational taxonomic units (OTUs) that were created based on a 97% 
sequence similarity cut-off value. 
doi:1 0.1 371 /journal.pone.0097629.g001 



highest OTU richness was observed for ITS86F/ITS4 with an 
average of 62 OTUs per sample (min = 42; max= 106). ITS IF/ 
ITS2 yielded on average 32 OTUs per sample (min=15; 
max = 60), whereas ITS3/ITS4 resulted in an average of 50 
OTUs per sample (min = 27; max =76) (Fig. 2B). Diversity was 
estimated with the inverse Simpson index. The inverse Simpson 
index differed significantly between ITS86F/ITS4 and ITS IF/ 
ITS2, whereas with ITS1F/ITS2 a lower diversity was found than 
with ITS86F/ITS4 (p = 0.04). However, no significant differences 
were found between ITS3/ITS4 and ITS1F/ITS2 or between 
ITS3/ITS4 and ITS86F/ITS4 (p = 0.31 and p = 0.53, respective- 
ly) (Fig. 2C). The average Good's coverage per sample obtained 
for ITS1F/ITS2 was 96.8% (min = 93.8%, max = 98.9%), where- 
as the average Good's coverage obtained for ITS3/ITS4 and 
ITS86F/ITS4 was 96.5% (min = 93.2%, max = 99.0%) and 
97.5% (min = 95.3%, max = 99.6%) respectively (Fig. 2D). Signif- 
icant differences in Good's coverage were found between ITS3/ 
ITS4 and ITS86F/ITS4 (p<0.01). However, no significant 
differences were found between ITS1F/ITS2 and ITS3/ITS4 
(p = 0.81) or between ITS1F/ITS2 and ITS86F/ITS4 (p = 0.31). 

Community similarity compared between primer pairs 

To compare the fungal community characterized with ITS 1 F/ 
ITS2, ITS3/ITS4 and ITS86F/ITS4 at the species and phylum 
level, a representative sequence of each OTU (as selected by 
Mothur) was blasted against the UNITE and INSD databases 
using the massBLASTer tool in PlutoF v2.0 [31]. Relative 
frequency distributions of the obtained species-level OTUs and 
phyla were analysed with chi-squared tests for the different primer 
pairs, based on the average abundances across replicates (four) and 
samples (seven). Representative reads of OTUs that could not be 
coupled to an accession of either the UNITE or INSD databases 
were considered as unidentified OTUs (indicated as not applicable 
"NA" in Table S3). A total of 5 1 unidentified OTUs were found of 
which 50 were found with ITS86F/ITS4 and 1 with ITS3/ITS4. 
BLAST scores and corresponding E-values for all OTUs can be 
found in Table S3. At the species level, differences were observed 
between the fungal communities identified by the three primer 
pairs studied (p<0.01). To give an idea of the fungal communities 
identified by each primer pair, pie charts displaying the top ten 
most abundant OTUs were constructed covering 68%, 62% and 
64% of all sequences obtained with ITS1F/ITS2, ITS3/ITS4 and 
ITS86F/ITS4, respectively (Fig. 3). Using the ITS1F/ITS2 primer 
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Figure 2. Parametrical comparison between the three primer pairs used in this study (ITS1F/ITS2, ITS3/ITS4 and ITS86F/ITS4). A. 

Average number of sequences obtained after quality trimming. B. Average number of operational taxonomic units (OTUs), based on a 97% sequence 
similarity cut-off value. C. Average inverse Simpson index. D. Average Good's coverage. Averages were calculated across replicates (four) and samples 
(seven) for each primer pair. Differences at the 95% significance level are indicated with an asterisk 
doi:1 0.1 371 /journal.pone.0097629.g002 



pair (targeting the ITS1 region) a total of 183 OTUs across all 
samples were observed, with the most abundant OTUs corre- 
sponding to Sistotrema sp. Fr. (27%), Rhizopogon luteolus Fr. (9%), 
Wilcoxina mikolae (Chin S. Yang & H.E. Wilcox) Chin S. Yang & 
Korf (8%), Cladophialophora minutissima M.L. Davey & Currah (7%), 
and Capronia sp. Sacc. (5%) (Fig. 3A). The primer pairs ITS3/ITS4 
and ITS86F/ITS4 (targeting the ITS2 region) identified 333 and 
346 OTUs across all samples, respectively. In line with ITS IF/ 
ITS2, the fungal communities identified with ITS3/ITS4 and 
ITS86F/ITS4 were also dominated by Sistotrema sp. (21-19%), but 
the subdominant OTUs were not exactly the same (Fig. 3B,C). 
Interesting to note is that the soil samples were dominated by 
ectomycorrhizal and ericoid mycorrhizal fungi and mycobionts 
from lichens. Based on field observations, we assumed that the 
fungal community in the pioneer forest that was sampled in this 
study would be relatively species poor compared to old forest soils 
[5] and that biotrophic fungi would dominate over saprotrophic 
ones. These assumptions were confirmed by all three primer pairs 
(Fig. 3). At the phylum level, differences in community compo- 
sition were found between all primer pairs tested (p<0.01 for all 
comparisons) (Fig. 4). Nevertheless, the majority of OTUs 
identified by all tested primer pairs belonged to the phyla 
Ascomycota (56% to 71%), followed by Basidiomycota (14% to 
17%). A minority of OTUs identified, belonged to the Zygomy- 
cota (3% to 4%), Chytridiomycota (3% to 4%) and Glomeromy- 
cota (0% to 3%) (Fig. 4). 

Repeatability of metabarcoding results 

The repeatability of the molecular identification of fungal 
OTUs from environmental samples was compared between the 



three tested primer pairs to assess their experimental robustness. 
Replicates of samples were compared for each primer pair using 
NMDS with Bray-Curtis dissimilarities. In this analysis, samples 
with a similar OTU-composition will have smaller Bray-Curtis 
distances than samples with more dissimilar OTU compositions. 
In general, for all three primer pairs, replicates from the same 
sample grouped closely together (especially for ITS3/ITS4) (Fig. 
SI). Hence, the results of molecular identification of fungal OTUs 
are fairly consistent between replicated samples using the current 
experimental set-up. In order to test the possibility that some 
OTUs are missed in metabarcoding analyses based on the 
amplification and sequencing of target DNA from a single DNA 
extraction, results from the four replicated DNA extractions of the 
same sample were compared (Fig. S2). This assessment was 
performed for the four most abundant OTUs, representing 
Sistotrema sp., Rhizopogon luteolus, Cladophialophora minutissima and 
Wilcoxina mikolae. From Fig. S2, it is clear that in some replicated 
extractions of the same sample abundant OTUs can be missed. 
These results indicate that PCR amplification and sequencing can 
best be performed on multiple DNA extractions from the same 
environmental sample that are pooled prior to PCR in order to 
obtain an accurate picture of a fungal community. 

Efficiency of primer pairs studied 

To test the amplification efficiency of the three primer pairs in a 
heterogeneous pool of DNA (environmental sample) a qPCR 
experiment was conducted. More specifically, a 2-fold dilution 
series, ranging from 1:1 to 1:64 dilutions of twelve randomly 
selected DNA samples were amplified with randomly selected 
ITS1F/ITS2, ITS3/ITS4 and ITS86F/ITS4 primers with MIDs 
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■ Sistotrema sp (ECM) 

■ Rhizopogon luteolus (ECM) 

□ Wilcoxina mikolae (ECM) 

■ Cladophialophora minutissima (END) 

■ Capronia sp (SAP) 

□ Chaenothecopsis pusiola (LICH) 

□ Epacris microphylla root assoc. fungus (ERM) 

□ Thelephora terrestris (ECM) 

■ Rhizoscyphus ericae (ERM) 

□ Lachnum sp (SAP) 

□ Remaining taxa 

■ Sistotrema sp (ECM) 

□ Cladonia fimbriata (LICH) 

□ Emericella nidulans (SAP) 

□ Wilcoxina mikolae (ECM) 

■ Cryptococcus podzolicus (SAP) 

■ Cladophialophora minutissima (END) 

■ Rhizoscyphus ericae (ERM) 

□ Oidiodendron reticulatum (SAP) 

□ Tricladium chaetocladium (SAP) 

■ Rhizopogon luteolus (ECM) 

□ Remaining taxa 

■ Sistotrema sp (ECM) 

■ Cladophialophora minutissima (END) 

□ Microsphaeropsis arundinis (SAP) 

□ Oidiodendron tenuissimum (SAP) 

□ Wilcoxina mikolae (ECM) 

■ Capronia sp (SAP) 

■ Rhizopogon luteolus (ECM) 

■ Cryptococcus podzolicus (SAP) 

n Cadophora finlandica (ECM/ERM) 

■ Rhizoscyphus ericae (ERM) 

□ Remaining taxa 



Figure 3. Relative abundance for the top ten most abundant species-level operational taxonomic units (OTUs), based on a 97% 
sequence similarity cut-off value, obtained for each of the three primer pairs studied (ITS1 F/ITS2, ITS3/ITS4 and ITS86F/ITS4). Reads 
that did not result in a BLAST hit against the UNITE or INSD databases were indicated as "not applicable (NA)". Ecological functions of OTUs are 
indicated between brackets behind the OTU identities (ECM: ectomycorrhizal, ERM: ericoid mycorrhizal, SAP: saprotrophic, LICH: lichenized, END: 
endophytic). OTUs not belonging to the top ten most abundant OTUs were pooled in the category "Remaining taxa". OTUs that appear exclusively in 
a single chart are indicated in grayscale. OTUs that can be found in multiple pie charts are indicated in colour. OTU abundance scores were averaged 
across replicates (four) and samples (seven). A. ITS1F/ITS2. B. ITS3/ITS4. C. ITS86F/ITS4. 
doi:1 0.1 371 /journal.pone.0097629.g003 



and 454 adaptors attached. For ITS1F/ITS2, exponential 
amplification was obtained between 24 and 32 PCR cycles for 
ten out of twelve samples (data not shown). For two samples no 
exponential amplification phase was obtained within 40 cycles 
with this primer pair. ITS3/ITS4 showed exponential amplifica- 
tion after 22 to 36 cycles for all twelve samples, whereas ITS86F/ 



ITS4 already showed an exponential amplification phase after 20 
to 31 cycles for all samples (data not shown). Average PCR 
efficiencies (± standard error) were calculated to be 76% (±4%) 
for ITS3/ITS4, 82% (±5%) for ITS1F/ITS2 and 97% (±6%) for 
ITS86F/ITS4 (Table 2). 
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■ Ascomycota 

■ Basidiomycota 

■ Chytridiomycota 

□ Glomeromycota 

■ Zygomycota 

□ NA 



■ Ascomycota 

■ Basidiomycota 

■ Chytridiomycota 
n Glomeromycota 

■ Zygomycota 
□ NA 



■ Ascomycota 

■ Basidiomycota 

■ Chytridiomycota 

□ Glomeromycota 

■ Zygomycota 

□ NA 



Figure 4. Relative number of OTUs belonging to different 
fungal phyla. OTUs that could not be assigned to a phylum 
were grouped together under "not applicable (NA)". Averages 
were calculated across replicates (four) and samples (seven). A. ITS1F/ 
ITS2. B. ITS3/ITS4. C. ITS86F/ITS4. 
doi:1 0.1 371 /journal.pone.0097629.g004 



Phylum-level PCR bias 

qPCR amplification efficiency did not significantly differ 
between primer pairs tested (ITS1F/ITS2, ITS3/ITS4 and 
ITS86F/ITS4), nor between phyla (Ascomycota, Basidiomycota, 
Glomeromycota and Zygomycota) (Fig. 5). Two-way ANOVA 
resulted in p = 0.14 for phylum and p = 0.59 for primer pair. 
Primer pair - rDNA target combinations with poor PrimerPros- 
pector scores tended to have slightly lower PCR efficiencies, but 
these differences were not significant. Strains of species used for 
this experiment and PCR efficiencies can be found in Table S4 
and Table S5, respectively. 

In silico evaluation of primers 

In a final analysis, the primer-to-target mismatches of the three 
primer pairs used in this study were evaluated with PrimerPros- 
pector [34]. PrimerProspector was used to calculate a score for 
each primer based on mismatches between primers and target 
DNA sequences. The closer the score of a primer is to 0, the fewer 
mismatches were detected between primers and target sequences. 
The average scores (± standard error) for primers used in our 
study were: ITS1F = 4.55 (±0.05), ITS2 = 0.70 (±0.03), 
ITS3 = 0.58 (±0.03), ITS4 = 3.96 (±0.04) and ITS86F = 0.52 
(±0.02) (Table 3). Moreover, it was found that 44% of the tested 
sequences had a mismatch with the last base at the 3' end of 
primer ITS IF. This particular mismatch between the last base at 
the 3' end of a primer sequence and a target sequence occurred 
with only 9%, 4%, 16% and 3% of the tested sequences for ITS2, 
ITS3, ITS4 and ITS86F respectively (Table 3). For comparison, 
also the primers suggested by Ihrmark et al. [35] and Toju et al. 
[36] were tested with PrimerProspector. Also in this analysis, 
ITS86F was found to have the best primer score of all tested 
primers (Table S6). 

Discussion 

Amplification and sequencing of short, standard DNA regions 
(metabarcoding) is becoming an increasingly popular tool for the 
characterization of fungal communities. Nevertheless, in most 
fungal metabarcoding studies, primers are generally used without 
being tested for their efficiency to amplify heterogeneous DNA 
pools, which may affect our view on studied fungal communities. 
Whereas the most commonly used primers in fungal metabarcod- 
ing studies were designed in the 90 s for species identification of a 
limited number of focal species, environmental metabarcoding 
studies generally aim to characterize diverse communities in 
environmental samples. Hence, primers used for fungal metabar- 
coding should be able to amplify a broad range of target DNA 
sequences in a sample that is also rich in non-target DNA and that 
may contain environmental contaminants [39]. Even though 
recent efforts have resulted in new primers that could amplify a 
large proportion of target fungal DNA sequences [35,36], an 



Table 2. Average PCR 


amplification 


efficiencies obtained for twelve environmental DNA samples using quantitative real-time PCR. 




Primer pair 


ITS1FJTS2 


ITS3JTS4 


ITS86FJTS4 


Average (%) 


82 


76 


97 


Standard error (%} 


4 


5 


6 


Minimum (%) 


64 


67 


78 


Maximum (%) 


97 


103 


120 


doi:1 0.1 371 /journal.pone.0097629.t002 
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Figure 5. Phylum-level PCR bias assessed using qPCR. Average PCR efficiencies were calculated for each phylum using 5 basidiomycetes, 5 
ascomycetes, 2 glomeromycetes and 3 zygomycetes. Error bars represent standard errors. No significant differences between primer pairs and phyla 
were found at the 95% significance level. 
doi:1 0.1 371 /journal.pone.0097629.g005 



experimental evaluation of PCR efficiency and primer perfor- 
mance should be performed on real environmental samples. 

Initially, also ITS1F/ITS86R was included in our study design, 
but this primer pair was discarded from the study as no 
amplification was obtained in exploratory PCR and gel-electro- 
phoresis tests. A plausible explanation for this failure can be found 
in the fact that the reverse primer (ITS86R) used and reported by 
Turenne et al. [27] and Vancov and Keen [26] contains an 
incorrect base at the 3' end of the primer sequence. In order to be 
the perfect reverse complement of ITS86F, the sequence of 
ITS86R should be 5 ' -TTC AAAG ATTC G ATGATTC AC - 3 ' , 
and not 5 ' -TTC AAAG ATTC GAT G ATTC AG- 3 ' as reported. 
GS-FLX sequencing of the amplicon pool resulted in 151,650 raw 
reads prior to quality trimming. Of these reads, 65,133 were 
assigned to their respective sample and 86,517 reads remained 
unassigned. The unassigned reads were investigated manually 
revealing that the majority were primer sequences probably 
resulting from primer dimers in our sequenced amplicon pool. 
Most likely, these primer dimers were not sufficiendy removed 
during post-PCR clean-up steps. 

Rarefaction curves were constructed for each primer pair 
(Fig. 1). These rarefaction curves indicate that the highest rarefied 
OTU richness and best coverage of the fungal community can be 
expected for the ITS86F/ITS4 and ITS3/ITS4 primer pairs. The 
average observed number of reads and the average observed 
number of OTUs (derived from these reads at a 97% sequence 
identity cut-off) indeed were highest for the ITS86F/ITS4 primer 
pair (797 reads and 62 OTUs on average per sample) and the 
ITS3/ITS4 primer pair (523 reads, 50 OTUs) and were much 
lower for the ITS1F/ITS2 primer pair (356 reads and 32 OTUs) 
(Fig. 2). The average observed diversity per sample, estimated by 



the inverse Simpson index, did not differ between ITS3/ITS4 and 
ITS86F/ITS4, but was significantly lower for ITS1F/ITS2 (Fig. 2). 
Overall, the low number of OTUs per sample found in the current 
study, are in sharp contrast with the more than 1000 OTUs per 
gram of forest soil found by Buee et al. [5] based on amplification 
with the ITS1F/ITS2 primer pair. This difference in richness may 
be explained by the fact that pioneer forests probably contain 
relatively fewer fungal species compared to old forest soils [5]. 
Additionally, overestimation or underestimation of species richness 
can also originate from data handling and analysis [40]. Based on 
the in silico performance and high Good's coverage calculated for 
ITS86F, it can be expected that the 62 OTUs found on average 
per sample by the ITS86F/ITS4 primer pair is close to the real 
species richness for the pioneer ecosystem growing on stabilised 
sand dunes which was studied here. The 50 OTUs per sample 
found by ITS3/ITS4 and the 32 OTUs found by ITS1F/ITS2, 
are probably underestimations due to a more narrow primer 
spectrum and/or lower PCR efficiencies. The fact that a high 
Good's coverage was found for the ITS1F/ITS2 primer pair 
despite a low observed OTU richness indicates that this primer 
pair is unable to multiply the ITS 1 region of a large number of 
fungi. This is also supported by the in silico analysis. In this analysis, 
ITS IF was shown to have the poorest primer score of 4.6 and its 
sequence was shown to have a mismatch at the final base at the 3' 
end of the primer (having a detrimental effect on amplification 
efficiency [41]) with no less than 44% of the tested fungal 
sequences (Table 3). The large number of mismatches between the 
ITS IF primer and its target sequences was previously also 
addressed by Bellemain et al. [42] and Ihrmark et al. [35]. In 
comparison, the ITS4 primer was given a score of 4.0 and was 
found to have a primer-to-target mismatch at the 3' end of the 



Table 3. Results of in silico testing of primers using PrimerProspector 1.0.1 [34]. 





Primer 


Number of sequences tested 


3' end base mismatch (%) 


Average score ± SE 


ITS IF (F) 


3748 


44% 


4.6 ±0.05 


ITS2 (R) 


4421 


9% 


0.7±0.03 


ITS3 (F) 


4421 


4% 


0.6±0.03 


ITS4 (R) 


4270 


16% 


4.0±0.04 


ITS86F (F) 


4421 


3% 


0.0±0.00 



Primers are indicated as forward (F) or reverse (R). Average PrimerProspector scores are shown ± standard errors (SE). 
doi:1 0.1 371 /joumal.pone.0097629.t003 
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primer with only 16% of the tested sequences. For the ITS2, ITS3 
and ITS86F primers a score of 0.7, 0.6 and 0.0 was obtained 
respectively (Table 3). These primers were shown to have a 
mismatch at the 3' end of the primer with only 9%, 4% and 4% of 
the tested sequences, respectively (Table 3), illustrating their broad 
amplification potential. Furthermore, our in silico analyses indicat- 
ed that the primers suggested by Ihrmark et al. [35] and Toju et al. 
[36] had more mismatches to their respective target sequences 
than ITS86F. 

To test how these parametrical differences would translate to 
amplification efficiency during PCR amplification preceding 
emulsion PCR (emPCR) and pyrosequencing, a first qPCR 
experiment was conducted. To this end, DNA was extracted 
from 12 soil samples and amplified with the same primer pairs 
used in the pyrosequencing experiment. The calculated PCR 
efficiencies were 82% for ITS1F/ITS2, 76% for ITS3/ITS4 and 
97% for ITS86F/ITS4 (Table 2). From these PCR efficiencies, it is 
clear that ITS86F/ITS4 amplified its target ITS regions with 
greater efficiency than the other two primer pairs. Contrary to our 
expectations from the in silico analysis, ITS3/ITS4 obtained a 
lower efficiency than the ITS1F/ITS2 primer pair. This could be 
explained by the fact that also other factors determine the 
amplification efficiency of PCR reactions beside binding and 
dissociation of primers to their target DNA sequences. Such 
factors include the temperature-dependent properties of target 
DNA sequences and primer sequences in the PCR mixture, the 
temperature-dependent behaviour of the used polymerase enzyme 
mixtures, the use of ROX as an endogenous reference dye, etc. 
[39]. Alternatively, the range of target sequences that ITS IF and 
ITS2 may bind to during PCR amplification is smaller, but the 
sequences that do get bound by these primers are amplified 
efficiendy. 

To see whether differences in amplification efficiency between 
primer pairs would also be reflected in the identities of the OTUs 
identified in the 454 amplicon pyrosequencing experiment, a 
representative read for each OTU was blasted against the UNITE 
and INSD databases and the BLAST hits with the highest score 
and a species or genus name were used to reconstruct the fungal 
community for each primer pair (Fig. 3). According to all three 
primer pairs, the soil fungal community was dominated by an 
OTU corresponding to Sistotrema sp. Additionally, all primer sets 
produced a number of OTUs that were commonly identified by all 
primer pairs (Fig. 3). The community identified by the three tested 
primer pairs still differed significantly, however. These differences 
confirm the finding that targeting either the ITS1 or the ITS2 
region may result in different pictures of the fungal communities at 
the OTU level, as was previously assessed by both in silico [43] and 
sequencing studies [40,44] . In addition, it was found that primers 
targeting the same ITS region do not necessarily result in the same 
OTU composition (Fig. 3), highlighting the importance of primer 
choice in a given study. However, it needs to be noted that in 
comparative studies, it has been shown that an ecological signal 
can be much stronger than the differences in community 
composition originating from primer choice [44]. 

At the phylum level, significant differences between ITS IF/ 
ITS2, ITS3/ITS4 and ITS86F/ITS4 were found as well (Fig. 4). 
Although in varying proportions, all three primer pairs identified 
more OTUs belonging to ascomycetes (70%, 71% and 56% 
respectively) than basidiomycetes (17%, 14% and 15%), but also 
Chytridiomycota (3%, 4%, 4%), Glomeromycota (0%, 3% and 
2%) and Zygomycota (3%, 3% and 4%) were detected (Fig. 4). 
This might suggest that more ascomycetes were present in the soil 
at the time of investigation. However, amplification of DNA from 
ascomycetes may be favoured relative to amplification of DNA 



from basidiomycetes as the ITS sequences for ascomycetes are 
generally shorter than basidiomycete ITS sequences (this is 
especially true for the ITS2 region [12]) and amplification of 
shorter DNA fragments is favoured during PCR. Whereas in 
previous in silico analyses indeed a phylum-level bias was expected 
for some of the primers used [42] , no such bias was found in the 
current study based on experimental data derived from qPCR of 
DNA extracted from 15 fungal species belonging to the major 
fungal phyla (Fig. 5). 

Whatever the aim of a metabarcoding study, results obtained 
from metabarcoding need to be reliable. To assess the repeatabil- 
ity of the fungal metabarcoding experiment, we analysed four 
replicate DNA extractions of seven soil samples separately. The 
analysis of all replicates of samples revealed that replicated analysis 
of the same sample with a specific primer pair generally results in 
similar fungal community compositions (Fig. SI). This is especially 
true for the ITS3/ITS4 and ITS86F/ITS4 primer pairs as their 
replicated samples clustered nicely together. However, this is less 
true for the ITS1F/ITS2 primer pair, where replicates of samples 
tend to have greater projected Bray-Curtis distances (Fig. SI). 
Moreover, we have shown that it is possible to miss certain OTUs, 
even abundant ones, when one sequences amplicon pools that are 
constructed from a single DNA extraction (Fig. S2). It is therefore 
advisable to extract DNA from environmental samples in multiple 
replicates, pool the eluates and perform PCR and sequencing on 
the DNA from the mixed eluate. This observation is in line with 
other studies performed previously, demonstrating that at least 
three replicated extractions are required to obtain a DNA pool 
that is representative for the microbial community present in a 
given soil sample [45,46]. 

Apart from the technical issues that were addressed in this study, 
our data also provided a glimpse at the fungal community present 
in the studied site. Based on field observations of above-ground 
basidiocarps, we assumed that pioneer pine forests in the Campine 
region in Belgium are dominated by biotrophic species (mosdy 
lichens, ectomycorrhizal and ericoid mycorrhizal fungi) over 
saprotrophic species. All three primer pairs confirmed this 
assumption, but they found different fungal OTUs to be 
dominant. According to the results obtained with ITS1F/ITS2, 
the fungal community in the studied site was dominated by OTUs 
corresponding to Sistotrema sp. (27%), followed by Rhizopogon luteolus 
(9%), Wilcoxina mikolae (8%) and Cladophialophora mimtissima (7%) 
(Fig. 3) [47]. These OTUs were also found to be very important 
members of the studied community according to ITS3/ITS4 and 
ITS86F/ITS4 as they appeared in the top ten of the most 
abundant OTUs found by both primer pairs, although in varying 
proportions (Fig. 3). Sistotrema sp., likely an important member of 
our studied ecosystem, was recently shown to be polyphyletic, 
containing both ectomycorrhizal and saprotrophic taxa [48] . The 
reads that were found in the current study correspond to Sistotrema 
strains that were sampled from ectomycorrhizal root tips of Pinus 
contorta Dougl. growing on coastal sand dunes [49]. This genus 
provides a fine example of the power of molecular tools, such as 
DNA metabarcoding, to draw attention to ecologically important, 
cryptic fungal species. Based on field observations alone 
(basidiocarps observations and root tip morphotying), we never 
expected this genus to be so abundant in this pioneer ecosystem. 

Concluding remarks 

In many fungal metabarcoding studies universal primers from 
previous phylogenetic or ecological studies are used without first 
performing an evaluation of their spectrum and performance for 
high-throughput sequencing, potentially resulting in a biased 
description of fungal communities. Whereas in silico PCR analyses 
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on sequences retrieved from sequence databases may suggest 
promising primers [35,36], we showed that an experimental set-up 
to evaluate their usefulness in practice provides complementary 
information on the actual performance of the primers for high- 
throughput sequencing of environmental samples. Indeed, here we 
demonstrated that the choice of primers has a significant impact 
on how fungal communities are translated into OTU communities 
and subsequent data analysis. As such, before setting up large-scale 
sequencing experiments, we recommend to first test a number of 
promising primer pairs, e.g. selected with in silico analyses, under 
real PCR conditions for a subset of the samples under 
investigation. In case an in-depth characterization of a fungal 
community is desired, the use of more than one primer pair is 
advisable. We also showed that quantitative real-time PCR, 
evaluating the efficiency of selected primer pairs, may help in 
selecting the most efficient primer pairs. After all, using primer 
pairs that are not very efficient in amplifying DNA from an 
environmental sample will undoubtedly result in a low number of 
reads, and consequently in biased community descriptions. 

In this study, the primer pair ITS86F/ITS4, which amplifies the 
ITS2 region of the fungal rRNA operon, was shown to be the most 
suitable primer pair for the characterization of fungal communities 
with metabarcoding. This primer pair not only resulted in superior 
amplification efficiency leading to a significandy higher number of 
reads, but also yielded a high number of OTUs belonging to 
different phyla. In addition, this primer pair resulted in a robust 
amplification reaction for the broadest range of samples and across 
replicated extractions. 
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Figure SI NMDS comparing community dissimilarities 
(based on Bray-Curtis distances) between replicates of 
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